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In plasmas at very low temperatures formation of neutral atoms is dominated by collisional 
three-body recombination, owing to the strong ~ y- 9 / 2 scaling of the corresponding recombination 
rate with the electron temperature T. While this law is well established at high temperatures, the 
unphysical divergence as T — > clearly suggest a breakdown in the low-temperature regime. Here, we 
present a combined molecular dynamics-Monte-Carlo study of electron-ion recombination over a wide 
range of temperatures and densities. Our results reproduce the known behavior of the recombination 
rate at high temperatures, but reveal significant deviations with decreasing temperature. We discuss 
the fate of the kinetic bottleneck and resolve the divergence-problem as the plasma enters the 
ultracold, strongly coupled domain. 



I. INTRODUCTION 

Since the first creation of ultracold plasmas (UCPs) 
via photo-ionization of laser-cooled atoms [JJ or cold 
molecules [2] , these systems prove to provide a well suited 
platform to study a range of plasma physics phenomena, 
such as collective waves 3 9, plasma expansion into vac- 
uum [TMT?1 , plasma instabilities [HI [TT5] and recombi- 
nation of neutral atoms [20H24] . Besides opening up a 
new parameter regime [251 - 12?] as well as promising ap- 
plications for nanotechnology [2"BH3"T] , UCPs offer a rare 
opportunity to study strongly correlated plasmas [32H35] 
in the laboratory. 

The degree of correlations is characterized by the 
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Coulomb coupling parameter V = a ^ BT , where e is 
the electron charge, fee is the Boltzmann constant, and 
a — (|7rp) _ s is the Wigner-Seitz radius for a plasma of 
density p. When the average potential energy ~ e 2 /a 
of the charges exceeds their thermal energy ~ k^T, i. 
e. when T > 1, the plasma is termed strongly cou- 
pled. Strong coupling phenomena studied in UCPs in- 
clude disorder-induced heating [35HJT] accompanied by 
kinetic energy oscillations [42-43] as well as liquid-like 
[33] and crystalline [45H45] behavior of the ionic plasma 
component. 

However, in contrast to dense strongly coupled plas- 
mas UCPs have a rather short life time on the order of 
several 10 2 /xs. The dominant decay mechanism is colli- 
sional recombination leading to the formation of highly 
excited neutral Rydberg atoms. The corresponding rate 
constant v ~ T~ 9 / 2 [31] has a strong dependence on the 
electron temperature T and can, hence, assume large val- 
ues in UCPs. In the underlying three-body recombina- 
tion (TBR) process two electrons collide in the vicinity 
of an ion to form a weakly bound Rydberg atom, where 
the remaining electron carries away the corresponding 
excess energy. Subsequent collisions of the formed Ryd- 
berg atom with free electrons can further de-excite the 
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atom, and the released binding energy leads to an in- 
crease of the free electron temperature T in the course of 
the plasma evolution HH H21 El ED EH ■ 

The strong temperature dependence of the recombi- 
nation rate can be readily established from simple scal- 
ing arguments. The collision frequency of an electron 
with average velocity v ~ \/T is given by v c = vb 2 p e , 
where p e is the electron density and b ~ T _1 is the char- 
acteristic distance of closest approach. Since the prob- 
ability of finding another electron within a distance b 
of the colliding pair is b 3 p e , one obtains a TBR rate 
of v ~ v c b 3 p e ~ p 2 T~z. In their seminal paper [45"] 
Mansbach and Keck presented a detailed study of TBR 
in ideal plasmas. Performing classical trajectory Monte- 
Carlo (CTMC) simulations of isolated three-body colli- 
sions and using rate equations they confirmed the T~^- 
scaling of the recombination rate and calculated the cor- 
responding proportionality constant. Subsequent studies 
based on CTMC calculations or rate equations 23, f5Tl - f53"] 
have since confirmed the strong temperature dependence, 
which was found to be in good agreement with experi- 
ments in hot and cold (T > 10000 K) plasmas j5Tl|521|54] 
as well as with measurements on UCPs with moderate 
coupling strength (r < 0.2) (22] [23]. 

On the other hand, the strong temperature dependence 
of the TBR rate ultimately suggests unphysically fast 
recombination in the ultracold regime. It was one of the 
major motivations of early UCP experiments [U [20] to 
shed light on this apparent divergence problem of the 
recombination rate. The conflicting timescales and the 
role of particle correlations become particularly evident 
by transforming to dimensionless units, employing the 
electronic plasma frequency ui p = \J kne 2 p e /m (m is the 
electron mass) such that 

v „fi T -9/2^ Wpr 9/2 (1) 

Consequently, for T > 1 the recombination rate v 
would become comparable or larger than the plasma 
frequency w p , whose inverse determines the typical 
timescale of electronic motion in the plasma. The un- 
physical crossover of timescales at T > 1 indicates a 
breakdown of the conventional recombination theory in 
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terms of three-body collisions and suggests a modifica- 
tion of the recombination rate in the strongly coupled 
regime. 

This fundamental problem has been addressed theoret- 
ically in a number of previous articles [53J I55H55] . In [551 - 
155] strong coupling in a dense plasma, where quantum ef- 
fects become important, has been considered and an en- 
hancement of the recombination rate coefficient (and ion- 
ization rate coefficient) due to correlation-induced contin- 
uum lowering has been found. On the other hand, studies 
of recombination in classical plasmas based on analytical 
estimates or numerical calculations [501 EH [Ml (Mj have 
found a suppression of recombination in the moderately 
to strongly coupled regime, but the proposed modifica- 
tions of the recombination rate yield different and con- 
tradictory results for the temperature scaling, such that 
the question of Rydberg atom formation in correlated 
plasmas still remains an unsettled issue. 

In this article, we investigate the classical recombi- 
nation of a single ion immersed in a strongly coupled 
electron plasma without employing any additional ap- 
proximation. Our calculations are based on Monte-Carlo 
(MC) sampling of classical molecular dynamics (MD) 
simulations that provide a natural extension of previ- 
ous three-body CTMC calculations [221 EFJ to account 
for strong electron-electron correlations and many-body 
interactions. Our results quantitatively reproduce the 
known behavior of the recombination rate [53] for T <C 1 
and are consistent with the results of Kuzmin and O'Neil 
[B"5] for the particular value of T = 0.6 studied in that 
work. We further discuss simulations of two-component 
plasmas for various initial configurations and temper- 
atures, that demonstrate strong disorder-induced elec- 
tron heating to T ~ 0.5. This conclusively excludes 
the possibility of a metastable plasma state with orders- 
of-magnitude suppression of recombination as suggested 
in [531 [M]. Nevertheless, strong coupling effects on the 
recombination rate are found to have observable conse- 
quences during the short-time evolution of UCPs. 

The article is organized as follows. Details of the 
plasma model as well as the numerical approach are given 
in section [nj where we also review the rate equation de- 
scription of TBR based on CTMC collision rates. In 
section [TIT] we discuss the obtained behavior of the bot- 
tleneck binding energy as function of T, allowing to cal- 



culate the recombination rate, presented in section IV 



Finally, section [V] provides simulation results for two- 
component neutral plasmas and a discussion of compet- 
ing heating effects and their consequences for the recom- 
bination dynamics. 



II. NUMERICAL APPROACHES 

In order to study Rydberg atom formation at a con- 
stant temperature and to isolate the recombination from 
other collision processes (see section [v]) we first consider 
the recombination dynamics of a single ion placed inside 



an electronic one component plasma (OCP), consisting 
of N electrons and a homogeneous neutralizing positively 
charged background. While this model can only provide 
a simplified description of a neutral two-component sys- 
tem (see section |v| it resembles the situation of antihy- 
drogen production experiments carried out at CERN [551 - 
155] . Here highly excited antihydrogen Rydberg atoms are 
formed via successive transits of antiprotons through an 
ultracold positron plasma. In these experiments atoms 
are formed predominantly via collisional recombination, 
which is however considerably modified by the presence 
of applied strong magnetic fields as has been extensively 
studied via CTMC calculations [B5H73) . 

Our simulations proceed in three steps. First, we equi- 
librate the electron OCP to a predefined temperature 
r _1 . Subsequently, we place a neutral atom at the cen- 
ter of the cubic simulation box, consisting of an ion and 
an electron at the origin with the electron having an ex- 
cess kinetic energy of ^ — . This procedure ensures that 
the total potential energy is not affected by the introduc- 
tion of the additional ion and at the same time gives a 
good description of atomic single-photon ionization used 
to produce UCPs pQ. Following the escape of the "photo- 
electron" we monitor the evolution of the surrounding 
plasma electrons. 



A. Molecular Dynamics Simulation 

The computationally most demanding part of the sim- 
ulation is the evaluation of the mutual interactions be- 
tween the N plasma electrons. Due the 0(N 2 ) of the 
corresponding numerical effort and the necessity to im- 
plement periodic boundary conditions (PBC) a straight- 
forward force calculation would be prohibitively demand- 
ing in view of the accuracy and ensemble size required for 
the present study. 

Both of these problems can be efficiently resolved by 
the fast multipole method (FMM) [75J which permits 
force calculations for large particle numbers with a com- 
plexity of O(N). The FMM algorithm divides the sim- 
ulation volume into a hierarchy of cubic subcells, and 
determines multipole expansions for the charge distribu- 
tion in each cell. The interaction of a certain particle with 
the particles in a distant cell can then be calculated much 
more efficiently and with a controllable error. Moreover, 
the FMM provides a natural implementation of PBC, as 
the periodic images of the simulation box can be treated 
as an upper extension of the hierarchy of cubic subcells 
[76] . The particular implementation of the FMM used in 
this work is described in [771 [75] . As discussed below, our 
algorithm also requires to calculate interactions among a 
small subgroup of particles. In this case, we perform a 
direct force summation and implement PBC via standard 
Ewald summation [751150] . 

In order to initiate the photo-electron at the central 
ion position we remove the singularity of the attractive 
electron-ion potential according to 



FIG. 1. (color online) Two dimensional projection of the cen- 
ter of the simulation cell with central ion (black dot) and the 
regions of different time steps. All electrons within the sphere 
r < Ro (red dots) are propagated with the same reduced time 
step At„ which is determined by the smallest electron - ion 
separation r m i n , according to R n < rwn < Rn-i (this ex- 
ample: At n — At 5). Electrons with r > Ro (blue dots) are 
propagated with the global time step Aip. 



V ion (r) = 



2r c y 1 3rfJ 



r >r c 
r < r r 



(2) 



where the softcore radius r c = 10~ 2 a was chosen suffi- 
ciently small to have no influence on the simulation re- 
sults, as has been checked by varying the value of r c . 

Consequently, the magnitude of the electron-ion force 
is limited by i^ ax = l/r 2 c . Typically, Fg ax will be 
much larger than the characteristic force between the 
electrons F ee ~ T -2 , which induces two vastly disparate 
timescales of the electron motion: Far away from the 
central ion electrons move comparatively slowly on a 
timescale ~ w p , while close to the ion scattered and 
bound electrons undergo a much faster dynamics. An 
efficient symplectic propagation of the electrons that ex- 
ploits different system time scales can be realized by 
the so-called reversible reference system propagator al- 
gorithm (r-RESPA) [EH [H2- In the present case, the 
principle idea is to evaluate the dynamics of distant elec- 
trons (blue dots in Fig. [I]) with a fixed, coarse-grained 
time step, while the electron motion in the immediate ion 
vicinity is resolved with a smaller, dynamically adapted 
time step. This guarantees an accurate description of the 
recombination process while maintaining computational 
costs at a minimum. 

Usually the r-RESPA split is either based on forces or 
on particles (52]. The two time scales discussed above 
call for a force based r-RESPA split by separating terms 
involving electron - ion interactions from terms which in- 
volve electron - electron interactions. However, as the 
condition F,^ ax ^> F ee holds only for those few elec- 
trons which are close to the ion, it is more efficient to 
additionally split the Hamiltonian based on particles, 
by separating terms involving electrons with positions 



FIG. 2. (color online) Schematics of the time step adaption 
as described in the text. 



r, < Ro = a from terms involving more distant electrons 
with ri > Rq. The total Hamiltonian H, describing the 
dynamics of the N electrons with position and mo- 
mentum Pi, (i = 1 . . . N), is thus split according to 



R = K f + K S + V f + V s + E 
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denote the kinetic energy and interactions of particles 
in the region of fast dynamics (r < Rq). Likewise, the 
remaining energies in the region of slow dynamics (r > 
Ro) are given by 



N n 2 
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(4) 



L i=0 



and Ebg is the particle - background interaction energy. 
The sum ^ L runs over all possible lattice vectors in the 
periodic lattice of image simulation boxes. Because of 
TBR collisions, which will take place in the region of 
fast dynamics, the electron - electron interaction among 
electrons in this region can become comparable to the 
ion - electron interaction and has thus been included in 

vf. 
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The majority of electrons participate in the slow dy- 
namics and are propagated with a fixed global time step 
Ato- The corresponding electron-electron forces includ- 
ing image charges are calculated within the FMM, while 
the electron-ion interaction is obtained by direct force 
summation combined with the Ewald summation to im- 
plement the PBC. Note that the forces arising from the 
image charges change only slowly and, hence, can be en- 
tirely accounted for in the slow dynamics (cf. eq. [4]) 

Each electron within the critical distance Rq = a par- 
ticipates in the fast dynamics. For appropriate time step 
adaption we divide this spherical region into concentric 
shells of radii R n = 2- n / 2 R , n = 1,2,... (see Fig. [l}. 
Once the smallest electron-ion distance falls below R n 
the time step for all electrons within R is decreased to 
At n+ i = 2- (n+1) At Q (see Fig. [l). For the electron prop- 
agation in the region of fast dynamics with the small 
time step At n the interaction Vt is calculated via direct 
force summation and the Ewald potential 7!) 80 K^llSoj . 
In order to avoid periodic changes of the time step At n , 
the time step is only decreased as long as there are elec- 
trons in the region r < R . Only when an electron leaves 
the region of fast dynamics its propagation is switched 
back to the global time step At (see Fig. |2j. We set 



At = 10" 



ensuring accurate propagation for both 



free and bound electrons. In fact, we achieve a very small 
maximum relative energy error in all simulations below 

0. 01%. We find that such a high accuracy is needed in 
order to obtain reliable converged results for the recom- 
bination rate. 

Formally our propagator for the fast electrons can be 
written as 

S* (Atn) = Uyf(^) U Kf (At n ) Uyf(^), (5) 

where we have used the notation employing propaga- 
tors UA(t) = e tVA of differential operators T>a — 
{q, A}. Here, {...} denotes the Poisson bracket, q = 
(ri, r2 . . . tn; Pi, P2 ■ • ■ Pn), and A stands for any of the 
energy terms of the system Hamiltonian The total 
system propagator is then given by 

S(Ato) = U V s(^)U K s(At ) (Sf(At n )f Uv«(^). 

(6) 

Besides the numerical challenges, a proper analysis of 
the MD data and in particular the identification and 
characterization of the formed classical bound states 
poses additional difficulties. The most straightforward 
way would be to monitor the electrons total energy E^, 

1. e. the sum of the i-th electron's kinetic energy and 
potential energy due to the ion, the remaining N — 1 
electrons, their periodic images and the positive homoge- 
neous background, and declare an electron bound when 
Ei < 0. While this criterion makes sense for isolated 
atoms, in the present case many-body interactions and 
possibly strong electron-electron correlations lead to a 
lowered and fluctuating ionization threshold [86] . In fact, 
already for moderate coupling strength one finds that 
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FIG. 3. (color online) (a) Typical time evolution of the min- 
imum electron energy i? m i n . Red parts correspond to bound 
electrons that are subsequently ionized, blue parts correspond 
to bound electrons that reach -E S mk without ionization for 
<j> c — 4-7T. (b) Corresponding energy densities pi on and p rec 
of the example trajectory (lines) and as obtained from the 
ensemble average (shaded). 



Ei < for several electrons. Those electrons can be con- 
sidered as weakly localized in slowly fluctuating potential 
wells formed by the correlated charges [S7]> but they are 
not bound by the ionic potential. Hence, once an electron 
becomes bound to the ion, its energy should drop signif- 
icantly below the energies of free electrons. Therefore, 
we use the lowest electron energy E m i n (t) = min(.Ej(i)) 

as an initial criterion to select a possibly bound electron, 
whose index is denoted by j m i n {t). 

A typical example of such an minimal energy trajectory 
Emin{t) is shown in Fig. [3^,. As can be seen, E m { n stays 
negative most of the time, and therefore provides only 
a necessary but not sufficient criterion for bound state 
assignment. In order to detect stable electron-ion orbits 
we adopt the procedure proposed in |88j and integrate the 
rotation angle <p of the j m i n th electron around the ion. If 
jmin changes its value, e.g. due to an exchange collision, 
integration starts again at zero. If the maximum rotation 
angle <fi m ax(jmin) of the j m inth electron exceeds a critical 
angle <p c the electron is considered to be bound. The 
binding energy E^(t) of a Rydberg atom corresponds to 
the parts in E m i n (t) where j m i n (t) fulfills this angular 
criterion: 



E h (t) 




1 'nit! x OminW) < <Pc 
(jmin(i)) > <t>c, 



(7) 



i.e. Eb(t) — in the absence of a bound state, as marked 
by the black segments in Fig. [3b. The influence of the 
precise value of <p c on the extracted recombination dy- 
namics will be discussed below. 

Following the initial capture, subsequent electron-atom 
collisions slightly (de)-excite the formed atom, lead to rc- 
ionization or occasionally drive the atom to significantly 
deeper binding energies. Such close collisions typically 
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cause electron exchange R)9"] and are marked by sharp 
peaks in E^{t) (see Fig. [3^,). 

The simulation is stopped as soon as E-^it) reaches a 
certain energy sink -E S ink- The energy sink was set to 
E sink /(e 2 /a) = -20 for r > 1 and to E sink /(e 2 /a) = -50 
for r < 1 to ensure that the probability for a re-ionization 
vanishes well before the energy sink is reached (see Fig. 
[4j and therefore does not affect the final results for the 
determined recombination dynamics. 

Physical quantities are extracted from averages over 
an ensemble of ~ 10 3 simulation runs per parameter set, 
produced in a Monte-Carlo sampling over the initial po- 
sitions and velocities of the electrons. 



B. Rate Equations 

To make direct comparison to previous CTMC calcu- 
lations, we also solved the corresponding rate equations 
for the recombination scenario discussed above. Here one 
calculates the evolution of level population densities p n 
of the recombining atom [STJ [5^1 GE3] according to: 



dt 



dt 



(t)R(n,n f )] 



+ p c (t) 3 R rcc (n) - p e (t) p„(t)R loa (n), I 
P C (t) [Pn'(t)R io n(n') - p c {t) 2 R, cc {n')\ 



(8b) 



We use the transition rates recently determined in [23] 
by CTMC calculations. The rate for excitation from the 
atomic level n to level n' is given by 



R{n, n) = k Q e 3 J, 2 e e ™ 
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+ 



9/2 



(e„ + 0.9) 7 /3 5/2 Ae4/3 



the rate for de-excitation from n to n' is given by 



(9) 



5/2 



R(n, n') — fco — 
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9/2 



(<*/+ 0.9)7/3 £ 5/2 Ae4/3 



, (10) 



the rate for recombination into an atomic level n is given 
by 



i? rec (n) 



11 ^/K/k B Tk n 2 A 3 

7/3 



e n +4.38ei- 72 + 1.32e n 
and the rate for ionization of atoms in level n is 
11 y/TZ/k B Tk e- £ " 



-Rion(n) 



7/3 



4.38 e}- 72 + 1.32e„ 



(11) 



(12) 



where k Q = e i / (k B TV^TZ) 1 e n = TZ/(n 2 k B T), Ae = 
|e n — e„'\, A = yj h 2 /(2tt m k B T) is the thermal de Broglie 
wavelength, 7Z is the Rydberg constant, and h is Planck's 
constant. 



III. KINETIC BOTTLENECK 

Generally, the recombination rate v is defined as the 
rate at which ground state atoms are populated in the 
plasma [3TJ |531 [521 US] • While such deeply bound states 
defy a classical description, it was shown in [35] that 
this rate can also be determined from the downward en- 
ergy flux through a kinetic bottleneck energy that di- 
vides weakly bound from stable atomic states. As the 
bottleneck typically lies in the classical region of binding 
energies this process can be described classically. 

The concept of the kinetic bottleneck is readily under- 
stood from the following simple arguments. Depending 
on its binding energy i?bn, a bound electron has a cer- 
tain probability Pi on {Eb) for collisional re-ionization and 
a probability to be successively driven to deeper binding 
energies until it eventually reaches the ground state with- 
out being re-ionized on its way. In our simulations the 
latter equals the probability P s i n k(-Eb) = 1 — Pion(Eb) for 
reaching the energy sink at E slnk , since Pion(-E'sink) ~ 0. 
The re-ionization probability decreases with deeper bind- 
ing and ultimately falls below the recombination proba- 
bility, such that the atomic states become more and more 
stable against ionizing electron-atom collisions. Hence, 
the kinetic bottleneck is defined as the energy £"b„ at 
which recombination starts to dominate, i.e. the energy 
at which 



-Pion(-Eb = Eb n ) = P s ink(^b — E\, 



(13) 



Three-body CTMC calculations predict a simple linear 
scaling of the bottleneck [4T)] 



Ehn 



-3.83 r 



(14) 



As the bottleneck energy is crucial for determining the 
recombination rate, we first need to check the validity 
of this simple law in the strong coupling regime. In the 
MD simulations, i?bn can also be determined from eq. 
( 13 1, where Pion{Eh) is calculated from the corresponding 



bound state energy densities obtained from the above 
described energy trajectories (see Fig. [3]) according to 



Ptot(e) 



6{E h (t)-e) dt 



(15) 



where r is the simulation time of a single simulation 
run and (...) denotes the average over the statistical 
ensemble. This total energy density can be split into 
two parts, Qt t(E h ) = g ion (E h ) + g ICC (E h ). g i0 n(Eb) con- 
tains only bound states that are subsequently ionized, 
while ftoc(^b) counts only energies of bound states that 
reach the energy sink without intermediate re-ionization 
(red and blue, respectively, in Fig. [3]d). The ioniza- 
tion probability Pion(-E'b) is then obtained from the ratio 
-Pion(-E'b) = g' t °°(^) 1 an( l shown in Fig. 4 for several val- 
ues of the Coulomb coupling parameterT. 
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FIG. 4. (color online) Ionization probability Pi on (£'b) for dif- 
ferent r. Symbols correspond to MD-data with <f) c = 4-7T. 
The lines show the ideal plasma prediction obtained by eq. 
(16 1. The crossing with P lon — | (dashed line) determines 
the location of the kinetic bottleneck £bn- 



For ideal plasmas (r — > 0) and within the adiabatic 
treatment of Bates, Kingston, and McWhirter [53], the 
ionization probability Pion(^b = —IZ/n 2 ) can be directly 
obtained from the collision rates eqs. (|9])-(12) [551 150"] 



Pion(^b) 



Rionjn) 

A(n) 



n' = l 



R{n,n') Rj on (n') 
A(n)A(n') 



■ (16) 



by summing over the probabilities of all possible path- 
ways in energy space that connect an atomic level of 
binding energy E\, — ~lZ/n 2 to the continuum, where 
A(n) = ^Sfl(n,n') + R ion (n) is the total rate for 
leaving level n. The first term in eq. ( 16 ) represents 



the probability that the bound electron will be ionized 
directly from level n. The second term accounts for an 
intermediate step via a level n' from which subsequent 
ionization occurs, and so on. 

The good agreement between our MD results and eq. 
(|16| in the regime of small to moderate T, shown in Fig. 
|4J attests to the accuracy of both approaches. In this 
regime the bottleneck energy can be straightforwardly 
determined according to eq. ( 13 ) and corresponds to the 
intersections of Pl on and the horizontal dashed line at 0.5 
in Fig. [4] At larger T-values, however, the ionization 
probability is suppressed to Pi on (Eb) < 0.5 over the en- 
tire range of binding energies, such that the bottleneck 
energy vanishes. 

The resulting temperature dependence of the bottle- 
neck energy is shown in Fig. [5]d. For small coupling 

parameters the MD simulations predict a linear scaling, 
_ 1 2 

i?bn = — 3.85T — , in quantitative agreement with the 
three-body CTMC result, eq. Q. However, for T > 2, 
the bottleneck drops to zero. In contrast to the ideal 
plasma case, where electrons still have to overcome the 
kinetic bottleneck barrier before recombination, stable 
atoms are formed directly in the strongly coupled regime. 
The disappearance of the kinetic bottleneck can be traced 
back to correlation-induced continuum lowering, which 



around T 2 leads to a merging of the ionization thresh- 
old and the bottleneck energy. To demonstrate this point, 
Fig. [5^ shows the fraction /< of free plasma electrons 
with a total energy of E.- L < — 4T _1 — . The simulation 
results yield a steep increase of /< around r ~ 2, at 
which the bottleneck, thus, has to disappear, in agree- 
ment with Fig. [5]d. Consequently, the critical angle C 
has almost no effect on the critical T at which the bot- 
tleneck energy drops to zero and can only slightly affect 
the value of E\, n for smaller T (see Fig. |5|t 



IV. RECOMBINATION RATE 

Having determined the location of the bottleneck we 
can now proceed to extract the recombination rate v from 
our MD simulations. This is done in a straightforward 
manner by calculating the time dependent recombination 
probability P rec (t), defined as the probability to observe 
a bound electron with binding energy E\>(t) < E^n at a 
time t. Fig. [6] shows examples of the obtained P rcc (t) for 
different coupling strength and a critical angle C = 47r. 
The numerical data is well fitted by an exponential bound 
state relaxation law of the form 



(17) 



which permits to extract the recombination rate v. 

Fig. [TJd shows the rate v as a function of the inverse 
critical angle <p c . One finds a linear dependence on 1/0 C! 
whose slope tends to increase with increasing coupling 
strength. The fact that in the considered range of <fi c 
all simulation results perfectly lie on a line, allows us to 
extrapolate to l/<f> c — > 0, corresponding to stable atomic 
states. 

For comparison with the weak coupling CTMC results, 




FIG. 5. (color online) (a) Fraction of free electrons with 
energy E < — 4T~ X — . (b) Bottleneck energy _Ebn as a func- 



tion of r for two different critical angles tf>, 
standard — 3.85F" 1 — scaling (black line). 
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FIG. 6. (color online) Recombination probability P lec {t) as a 
function of time for different coup ling strength T and <f) c = 4tt. 
The lines represent fits of eq. (17 1 to the MD simulation 
result. 



we also calculate the recombination rate 



v = [! - Pion(-n/n 2 )] R rcc (n) 



(18) 



as obtained from eq.(ll) and (16). Fig. [7^i shows the 
extrapolated many-body MD (circles) and three-body 
CMTC (squares) results for the recombination rate as 
a function of T. In the weak coupling regime we find 
good quantitative agreement with the T -9 / 2 scaling [23] 



v « 0.019 w p r 9/2 



(19) 



Notably, the MD results demonstrate the high accu- 
racy of the rate equation description even for moderate 
coupling strength V < 0.3, corresponding to typical pa- 
rameters in the long-time evolution of UCPs [TTJ 150] , 
In the regime of strong Coulomb coupling, however, 
one finds a significant suppression of the recombina- 
tion rate. Our result approaches a constant value of 
v ~ 0.03w p with increasing T, thereby resolving the ap- 
parent timescale paradox described in section [i] At in- 
termediate T-values, our results are consistent with pre- 
vious MD simulations of two-component plasmas |65j . 
These two-component simulations also predict a sup- 
pression by a factor of ~ 2 for the particular value of 
T = 0.6 studied in [65], suggesting that the present OCP 
model should provide a good description of recombina- 
tion in neutral plasmas. In this case, however, addi- 
tional disorder-induced electron heating [37] [38] due to 
the strong attractive electron-ion interaction limits the 
range of realizable Coulomb coupling parameters, as will 
be briefly discussed below. 



V. TWO-COMPONENT PLASMA 
SIMULATIONS 

We also performed MD simulations of a two component 
plasma with N ions and N electrons in a cubic simulation 



FIG. 7. (color online) (a) Recombination rate u as a function 
of inverse coupling strength T" 1 calculated with MD simula- 
tions (circles) compared to the F -9 ' 2 -scaling obtained by eq. 
( 18 1. The red line corresponds to eq. ( 19 1, the blue line serves 
as guide to the eye. (b) Dependence of v on inverse critical 
angle l/<j} c with fitted extrapolation (lines). 



cell with PBC. For these simulations, all interactions and 
the corresponding PBC are calculated within the FMM. 
We use a very small global time step At = 10~ 5 , ensur- 
ing an accurate treatment of even the lowest bound states 
observed in the simulations. In analogy to the previously 
described simulation scheme, the full plasma simulations 
start with N randomly distributed atoms which are pho- 
toionizcd at t = as detailed in section Hi Al The initial 
kinetic excess energy E — =S- — determines the initial 
effective coupling strength Tq, i.e. the scaled kinetic en- 
ergy r^ 1 . 

However, since this procedure creates a highly non- 
equilibrium plasma, it takes a finite time to establish a 
well defined electron temperature. In order to character- 
ize the corresponding initial relaxation we monitor the 
evolution of two different coupling parameters, defined 
through ith order momenta (v l ) of the free-electron ve- 
locity distribution 



r 2 



3/(v 2 ) 



(20a) 
(20b) 



In local equilibrium, i.e. once the electrons have es- 
tablished a Maxwellian velocity distribution, Ti = T2- 
Indeed the simulation results shown in Figs. [5] and [9] 
show that both definitions of T approach each other on 



a timescale 



such that one can speak of an electron 



temperature for t > lu~ 1 . 

On the same time scale, the electrons heat up due to 
disorder-induced heating [37l [38] . As shown in Figs. [8] 
and|9j for both very large (r = 50, Fig{8]) and moderate 
(r = 1, Figj9| initial effective coupling strength, the 
Coulomb coupling parameter relaxes to a value of T rj 0.5 
during the initial relaxation stage, which is a factor of ~ 2 
smaller than found in [35] but agrees with the findings of 
more recent MD simulations 1ST]. 
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FIG. 8. (color online) Evolution of the normalized free elec- 
tron density f e (a) and electronic coupling parameter V (b) 
calculated with MD simulation for <f) c — 4 n and rate equation 
(RE) for an initial excess energy corresponding to coupling 
strength of To = 50. 



FIG. 9. (color online) Same plot as in Fig. [8] but for To = 1. 



We have investigated the amount of initial heating for 
a range of initial conditions, including highly pre-ordered 
states, where ions and electrons have been placed on reg- 
ular lattice structures. In contrast to the ionic plasma 
component, where such a pre-ordering leads to signifi- 
cant suppression of the heating due to the repulsive ion- 
ion interactions [221 EES, the attractive electron- ion in- 
teraction is found to cause electron heating to T < 0.5 
irrespective of the initial state. This clearly excludes 
the existence of metastable, very strongly coupled two- 
component plasma states, in which recombination is sup- 
pressed by orders of magnitude, as has been suggested 
recently [63l [64] on the basis of numerical simulations. 

Nevertheless, the results of the previous sections have 
shown that the recombination dynamics differs signifi- 
cantly already for T s» 0.5, yielding a suppression of v by 
a factor of 2. A comparison to the time evolution of the 
free electron number obtained from the rate equations 
(see section II B ) reveals significant deviations, which in- 
crease with decreasing initial energy (see Figs. |9]and[8|. 
Moreover, despite the fact that T has relaxed to ~ 0.5 
after ~ 50 ui^ 1 for both To = 1 and To = 50, the number 
of weakly bound Rydberg atoms differs by about a factor 
of ~ 2, which we attribute to the longer relaxation time 
of bound states. We anticipate, that these effects and 
deviations from the traditional treatment of recombina- 
tion may be observable via short-time probing of UCP 
dynamics. 

Indeed, recent measurements of recombination fluores- 
cence on a sub-microsecond time scale suggest such devi- 
ations [24 . In this experiments, the time dependent fluo- 
rescence from low-lying transitions of recombined atoms 
has been measured with a time resolution < 100 ns. At 
high temperatures, i.e. in the weakly coupled regime, the 
initial signal S(t) was found to rise proportional to p^, 
consistent with the picture of isolated three-body colli- 
sions, for which S(t) ~ pvt oc pu} p tT 9 / 2 cx p\. For lower 



temperatures, detailed simulations based on three-body 



CTMC rates predict a density scaling 



while the 



experiment shows a scaling ~ p£ ■ Here it is interesting 
to note, that the recombination rate around r ~ I al- 
ready shows a different scaling v ~ WpF 3 (see Fig[7]), 
giving pvt oc p 2 / 5 . As described in [23], the fluores- 
cence signal is, however, determined also by the initial 
disorder-induced heating as well as heating due to the 
formation of Rydberg atoms themselves, such that this 
simple comparison should be regarded as qualitative only. 
Nevertheless, MD simulations, as described in this work, 
combined with a detailed treatment of the radiative cas- 
cade of deeply bound states to compare with such mea- 
surements may elucidate the role of correlation effects in 
recombination dynamics of UCPs. 




FIG. 10. (color online) Recombination rate v calculated with 
MD simulation (circles) for intermediate coupling strength T. 
At the onset of strong correlations the rate scales with F 3 
(line). 
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VI. SUMMARY 

We have presented extensive numerical simulations of 
Rydberg atom formation in plasmas, that take into ac- 
count correlations and many-body interactions between 
the plasma electrons. This allows to stretch the focus 
of such studies deep into the strongly coupled regime, 
beyond the range of validity of three-body CTMC calcu- 
lations. 

We find quantitative agreement for the recombination 
rate with previous rate equation calculations [23] in the 
weakly coupled regime. Such simplified treatments are 
shown to yield an excellent description, even for Coulomb 
coupling strengths of up to T — 0.3, which covers the 
typical T-values obtained in the long-time dynamics of 

ucpshiisd]. 

However, as the electron plasma becomes strongly cou- 
pled the bottleneck is found to disappear in the lowered 
continuum due to increasing electron-electron correla- 
tions. In this strongly coupled regime, v is shown to ap- 
proach a constant value well below the plasma frequency 
w p , resolving the temperature-divergence problem of the 
common three-body recombination rate, v ~ y-s/ 2 ) m 
the ultracold domain. 



MD simulations of two component plasmas show that 
the achievable coupling strength in neutral plasmas is 
limited to T ~ 0.5, in agreement with recent simula- 
tions discussed in while contradicting the findings of 
[SHIM]. Nevertheless, a comparison to the MD results for 
such coupling parameters suggest that deviations from 
common rate equation descriptions may be observable 
in the short-time dynamics of UCPs. We finally note, 
that recent experiments on molecular ultracold plasmas 
HZ] j realizing much higher densities than atomic sys- 
tems, show strong deviations from the expansion behav- 
ior of atomic systems |10j . which, thus far, has been well 
described within simple rate equation treatment of Ry- 
dberg atom formation [Til US EI] ■ Exploring the origin 
of these deviations, however, requires to account for ad- 
ditional molecular processes [M|, that may also alter the 
plasma expansion behavior. 
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